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I.  INTRODUCTION 


Experimental  testing  and  gasdynamic  modeling  of  the  ram  acceleration  technique  for 
in-bore  projectile  propulsion1 2 3 4'5  has  been  initiated  at  the  U.S.  Army  Ballistic  Research  Lab¬ 
oratory  (BRL)  under  the  Hybrid  Inbore  RAM  (HIRAM)  propulsion  program.6  Numerical 
solutions  of  the  Navier-Stokes  equations  have  been  obtained  via  computational  fluid  dy¬ 
namics  (CFD)  for  non-reacting  and  reacting  flows.7  Numerical  simulations  can  be  used  as 
a  system  design  aid  and  as  a  means  by  which  geometric  and  fluid  dynamic  scaling  phe¬ 
nomena  are  investigated.  Reacting  flows  considered  in  this  report  include  the  thermally 
choked,  subsonic  combustion  ram  acceleration  mode. 

For  non-reacting  flow  simulations  the  Rockwell  Science  Center  USA-PG  (Unified  Solu¬ 
tion  Algorithm  Perfect  Gas)  code8-9  is  used.  In  this  code  the  Navier-Stokes  equations  arc 
cast  in  conservation  form  and  converted  into  a  set  of  algebraic  equations  using  upwind  and 
central  finite  differences,  finite- volume  formulations  and  a  second-order  TVD  (total  vari¬ 
ation  diminishing)  scheme.  The  code  is  capable  of  simulating  mixed  subsonic/supersonic 
flowfields.  The  USA-PG  code  has  been  used  by  the  BRL  to  investigate  solid-fuel  ramjet 
projectiles.10  For  ram  accelerator  simulations  the  USA-PG  code  was  run  using  inviscid, 
laminar,  and  turbulent/separated  flow  modes  thus,  viscous  flow  effects  such  as  turbulent 
wall  boundary  layers  and  wake  flow  mixing  are  investigated.  Flowfield  shock  wave  gener¬ 
ation  and  reflection  patterns  are  computed.  Fluid  dynamic  scaling  phenomena  associated 
with  scaled-up  projectile  geometries  (i.e.  38mm  to  120mm)  are  investigated  as  well.  Com¬ 
puted  ram  accelerator  tube  wall  pressures  are  compared  to  data  measured  at  the  University 
of  Washington  using  a  38mm  system  and  inert  gas  mixtures. 

For  the  reacting  flow  simulations  the  BRL-LFD  RAMCOMB  (RAMjet  COMBustion) 
code11  (solid  fuel  diffusion-flame  technique)  was  modified  for  premixed  gaseous  fuel 


1  Bruckner,  A.P.,  Bogdanoff,  D.W.,  Knowlen,  C.,  and  Hertzberg,  A.,  “Investigation  of  Gasdynamic  Phenomena  Associated 
with  the  Ham  Accelerator  Concept,”  AIA  A-87- 1327,  Proceedings  of  the  AIAA  1 9th  Fluid  Dynamics,  Plasma  Dynamics  and 
Lasers  Conference,  June  8-10,  198 7,  Honolulu,  Hawaii, 

2  Hertzbcrg,  A.,  Bruckner,  A.P.,  and  Bogdanoff,  D.W. ,  “Ram  Accelerator:  A  New  Chemical  Method  for  Accelerating  Pro¬ 
jectiles  to  Lfltrahigh  Velocities.”  AIAA  Journal.  Vol.  26,  No.  2,  Feb.  1988,  pp.  195-208. 

3 Bruckner ,  A.P.,  Knowlcn,  C Scott,  K.A.,  and  Hertzberg,  A.,  “High  Velocity  Modes  of  the  Thermally  Choked  Ram 
Accelerator,”  A IA  A-88-2925,  Proceedings  for  the  AIAA  24th  Joint  Propulsion  Conference  July  11-13,  1968 ,  Boston,  MA. 

4  Burnham,  E.A.,  Kuil,  A.E.,  Knowlcn,  C.t  Bruckner,  A.P.,  and  Hertzberg,  A.,  “Operation  of  the  Ram  Accelerator  in  the 
Transdetonative  Velocity  Regime,”  AIAA-90-1985,  Proceedings  for  the  AIAA  26th  Joint  Propulsion  Conference ,  July  16-16, 
1990,  Orlando,  FL. 

' Bruckner ,  Hertzberg,  A.,  and  Knowlcn,  C.,  “ Review  of  Ram  Accelerator  Propulsion  Modes,”  Proceedings  of  the  27th 
JANNAF  Combustion  Subcommittee  Meeting,  November  5-9,  1990,  Warren  AFB,  Cheyenne,  WY. 

15  Kruczynski,  D.L.,  “Large  Caliber  Applications  of  Ram  Accelerator  Technology,”  Proceedings  of  the  27th  JANNAF  Com¬ 
bustion  Subcommittee  Meeting,  November  5-9,  1990,  Warren  AFB,  Cheyenne,  WY. 

7  Nusca.  M.J.,  “Numerical  Simulation  of  Reacting  Flow  in  a  Thermally  Choked  Ram  Accelerator,”  Proceedings  of  the  27th 
JANNAF  Combustion  Subcommittee  Meeting,  November  5-9,  1990,  Warren  AFB,  Cheyenne,  WY. 

8  Chakravarthy,  S.R.,  Szema,  K.Y.,  Goldberg,  U.C.,  Gorski,  J.J.,  and  Osher,  S.,  “ Application  of  a  New  Class  of  High 
Accuracy  TVD  Schemes  to  the  Navier-Stokes  Equations,”  AIAA-85-0165,  Proceedings  of  the  AIAA  23rd  Aerospace  Sciences 
Meeting,  Jan.  14-17,  1985,  Reno,  NV. 

9  Chakravarthy,  S.R.,  Szema,  K.Y.,  and  Haney,  J.W.,  “Unified  Nose-to-  Tail  Computational  Method  for  Hypersonic  Vehicle 
Applications,”  A  IA  A-89-2564,  Proceedings  of  the  AIAA  6th  Applied  Aerodynamics  Conference,  June  6-6,  1969,  Williamsburg, 
VA. 

10  Nusca,  M.J.,  Chakravarthy,  S.R.,  Goldberg,  U.C.,  “Computational  Fluid  Dynamics  Capability  for  the  Solid- Fuel  Ramjet 
Projecttle AIAA  Journal  of  Propulsion  and  Power,  Vol.  6,  No.  3,  May-June  1990,  pp.  256-262  (see  also  AlAA-87-2411 
and  BRL-TR-2958,  Dee.  1988). 

11  Nusca,  M.J.,  “Steady  Flow  Combustion  Model  for  Solid- Fuel  Ramjet  Projectiles,”  AIAA  Journal  of  Propulsion  and  Power, 

Vol.  6,  No.  3,  May-June  1990,  pp.  348-352  (see  also  AIAA-89-2797  and  BRL-TR-2987,  April  ! 969).  ' 
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combustion.  The  steady  2D/axisymmetric  Navier-Stokes  equations  are  written  in  stream 
function-vorticity  form  and  solved  using  a  Gauss-Seidel  relaxation  scheme.  These  equations 
include  conservation  of  chemical  species  and  reacting  flow  source  terms  in  the  energy 
equation.  A  one-step,  global,  temperature  and  mass  controlled,  finite-rate  reaction  between 
fuel  (CH4)  and  oxidizer  (O2)  is  considered.  The  products  of  this  reaction  (e.g.  CO2  and 
H2O)  are  lumped  into  one  specie  with  mass-averaged  thermodynamic  properties.  Diluent 
species  (reactants  and  products)  are  also  considered.  Reaction  rates  are  defined  using  the 
Law  of  Mass  Action  and  Arrhenius  expressions.  Constants  in  the  rate  expressions  as  well 
as  the  molar  specific  enthalpy,  heat  of  formation,  specific  heat,  and  molecular  weight  of  the 
species  are  extracted  from  previous  studies  of  hydrocarbon  reactions.  Upstream  subsonic 
flow  conditions  for  the  reacton  zone  are  obtained  using  the  USA-PG  code. 


II.  BACKGROUND 


The  ram  accelerator  technique  was  first  investigated  via  experimental  test  firings  at 
the  University  of  Washington  (UW).1'5  The  UW  ram  accelerator  facility  (Figure  1)  uses 
a  light  gas  gun  (e.g.  helium  driver  gas)  to  accelerate  projectiles  up  to  1300  m/s.  The 
muzzle  is  connected  to  a  perforated  tube  and  evacuation  tank  which  serve  as  a  dump  for 
the  driver  gas  prior  to  entrance  into  the  16m  ram  accelerator  tube.  This  tube,  which  can 
be  divided  into  sections  separated  by  diaphragms,  is  filled  with  a  pressurized  fuel/oxidizer 
mixture  and  is  instrumented  at  40  axial  locations.  Instrumentation  consists  of  pressure 
transducers,  fiber-optic  light  guides,  and  magnetic  transducers.  Thin  magnetic  sheets  are 
mounted  in  the  nose-body  joint  and  in  the  base  of  the  projectile  (see  Figure  1).  When 
the  projectile  passes  electromagnetic  transducers  on  the  accelerator  tube,  these  magnets 
induce  signals  that  are  used  to  determine  the  distance-time  history  (i.e.  velocity)  of  the 
projectile. 

The  projectile  consists  of  an  axisymmetric  cone-boattail  body  with  stabilizing  fins  to 
center  it  along  the  launch  tube  axis  (see  Figure  2,  fins  excluded).  The  projectile  is  shaped 
like  the  centerbody  of  a  ramjet  engine  and  is  injected  into  a  stationary  tube  filled  with  a 
pressurized  gas  mixture  of  hydrocarbon  fuel  (e.g.  CH4),  oxidizer  and  diluents  such  as  CO2, 
N2,  He,  and  Ar.  There  is  no  propellant  on  board  the  projectile.  The  tube  resembles  the 
outer  cowling  of  a  conventional  ramjet  engine.  When  the  injection  velocity  is  greater  than 
the  sound  speed  of  the  gas  an  oblique  shock  system  develops  on  the  projectile  which  can 
sustain  combustion  on  the  surface  and  aft  of  the  projectile.  In  this  way  the  energy  release 
process  travels  behind  the  projectile.  Thrust  is  generated  by  the  action  of  high  pressure 
reacting  gases  on  the  rear  part  of  the  projectile.  Various  combustion  ignition  mechanisms 
have  been  used.  The  initial  gas  pressure,  fuel/oxidizer  composition,  and  sound  speed 
can  be  selected  to  achieve  the  desired  acceleration  and  projectile  velocity  at  tube  exit. 
Diluents  are  used  to  tailor  the  acoustic  speed  of  the  mixture  so  that  the  initial  Mach 
number  of  the  projectile  exceeds  the  minimum  required  (~  2.8)  to  start  the  diffuser  (i.e. 
the  minimum  projectile/tube  clearance  at  maximum  projectile  body  diameter)  and  tailors 
the  heat  release  of  combustion  to  a  level  that  stabilizes  the  shock  system  on  the  projectile 
body.  Excessive  diluent  results  in  low  projectile  acceleration  levels,  whereas  insufficient, 
diluent  concentration  can  cause  pre-ignition  forward  of  the  projectile. 
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The  balance  of  forces  on  the  projectile  is  composed  of  the  drag  force  on  the  forebody 
and  a  thrust  force  produced  by  the  high  pressure  combustion  products  on  the  projectile 
afterbody.  Normally  a  net  thrust  is  obtained  since  the  pressure  of  the  combustion  products 
is  higher  than  the  compressed  gases  downstream  of  the  nose  shock.  Gas  mixtures  with 
higher  heat  of  reaction  yield  greater  net  thrust.  Combustion  must  begin  on  the  afterbody 
of  the  projectile  in  order  to  maximize  the  thrust.  If  combustion  starts  on  the  conical  nose, 
the  high  pressure  due  to  the  reaction  products  will  contribute  to  drag  forces  and  thus, 
will  reduce  the  net  thrust.  Similarly,  if  combustion  starts  aft  of  the  cone-boattail  junction 
on  the  projectile,  then  only  a  fraction  of  the  afterbody  area  will  be  exposed  to  the  high 
pressure  gases  and  the  thrust  will  be  below  its  full  potential  value. 

Several  modes  of  ram  accelerator  operation  have  been  investigated  including  thermally 
choked,  transdetonative,  oblique  detonation  wave,  and  superdetonative.  These  modes  are 
distinguished  by  operating  velocity  and  the  manner  in  which  the  combustion  process  is 
initiated  and  sustained.  The  thermally  choked,  subsonic  combustion  mode  will  theoret¬ 
ically  allow  projectile  acceleration  up  to  the  Chapman- Jouguet  (C-J)’  detonation  speed 
of  the  gaseous  propellant  mixture.  Operating  modes  that  involve  supersonic  combustion 
have  demonstrated  acceleration  beyond  the  C-J  speed.5  Numerical  simulation  of  ram 
accelerator  flowfields  has  been  documented  for  supersonic  flow  such  as  the  oblique  deto¬ 
nation  wave  and  superdetonative  modes.12'15  This  report  describes  CFD  modeling  of 
38mm  and  120mm  ram  accelerator  systems  for  non-reacting  and  reacting  gas  mixtures  in 
the  subsonic  flow  thermally  choked  mode. 

Demonstrations  of  thermally  choked  ram  acceleration  at  the  University  of  Washington 
have  achieved  velocities  up  to  2.5  km/s  in  a  38mm  diameter  tube  with  a  projectile  mass 
of  about  70g.  A  projectile  is  injected  into  the  accelerator  tube  at  700-1200  m/s.  The 
forebodv  angle  of  the  projectile  is  such  that  the  oblique  shock  system  does  not  initiate 
combustion.  A  normal  shock  is  located  downstream  of  the  maximum  cross-sectional  area 
of  the  projectile  (Fig.  3);  theoretically  this  shock  is  also  not  strong  enough  to  ignite  the  gas 
mixture3  provided  that  the  projectile  velocity  does  not  exceed  about  Mach  4.  Combustion 
is  initiated  immediately  behind  the  projectile  upon  entrance  to  the  ram  accelerator  tube 
(a  published  discussion  of  ignition  is  forthcoming  from  UW).  The  base  of  the  projectile  can 
act  as  a  flameholder.  The  combustion  zone  spreads  to  the  full  diameter  of  the  tube  and 
reaches  a  thermal  choking  condition  some  distance  downstream  of  the  projectile.  Thermal 
choking  and  the  heat  release  due  to  combustion  stabilizes  the  normal  shock.  In  principle 
the  concept  is  scalable  for  projectile  mass  in  kilograms.  Verification  of  the  scaling  principles 
is  underway  at  the  BRL  using  both  experiments6  and  numerical  simulations.' 

*  A  normal  shock  in  a  constant  area  duct  followed  by  heat  addition  and  thermal  choking  in  steady  flow  constitutes  a  C-J 
detonation  wave. 

12  Bogdanoff,  D.W.,  and  Brackett,  D.C.,  “ A  Computational  Fluid  Dynamics  Code  for  the  Investigation  of  Ramjet-In- Tube 
Concepts, ”  AIAA-87- 1978,  Proceedings  of  the  AIAA  23rd  Joint  Propulsion  Conference,  June  29-July  2,  1987,  San  Diego,  CA. 

Yungster,  S.,  Eberhardt,  S.,  and  Bruckner,  A.P.,  “Numerical  Simulation  of  Shock* Induced  Combustion  Generated  by  High * 
Speed  Projectiles  in  Detonable  Gas  Mixtures,”  AIA  A-89-0673,  Proceedings  of  the  AIAA  27th  Aerospace  Sciences  Meeting,  Jan. 
9*12,  1989,  Reno,  NV. 

14  Yungster,  S.t  and  Bruckner,  A.P.,  UA  Numerical  Study  of  the  Ram  Accelerator  Concept  in  the  Superdetonative  V 'elociiy 
Range,”  AI AA-89-2B77,  Proceedings  of  the  AIAA  25th  Joint  Propulsion  Conference,  July  10*12,  1989,  Monterey,  CA. 

15  Chuck,  C.,  and  Eberhardt,  S.,  “ Numerical  Simulation  of  Hypersonic  Oblique  Shock-  Wave/ Laminar  Boundary* Layer  In¬ 
teraction  with  Shock  Induced  Combustion,”  AIAA-90-0149,  Proceedings  of  the  AIAA  28th  Aerospace  Sciences  Meeting,  Jan. 
8-11,  1990,  Reno,  NV. 


III.  NON-REACTING  FLOW 


The  CFD  approach  described  in  this  section  can  be  used  to  predict  the  compress¬ 
ible  flowfield  around  aerodynamic  bodies  by  solving  two-dimensional  (2D)  and  three- 
dimensional  (3D)  Reynolds- averaged  Navier-Stokes  (RANS)  equations.  The  USA-PG2 
and  USA-PG3  codes  were  developed  by  Chakravarthy.8’9  The  Navier-Stokes  equations  are 
written  using  the  perfect  gas  assumption  however,  versions  of  the  code  for  a  real  gas  have 
been  documented  as  well.16’17  Both  laminar  and  turbulent  flows  can  be  investigated 
thus,  a  turbulence  model18  is  required  for  closure.  In  addition,  backflow  regions  can  be 
present  thus,  a  backflow  turbulence  model19  is  included.  The  equations  are  transformed 
into  conservation  law  form  and  discretized  using  finite- volume  approximations.  The  USA- 
PG  code  uses  a  class  of  numerical  algorithms  termed  total  variational  diminishing  (TVD) 
which  do  not  require  the  inclusion  of  explicit  smoothing  or  global  dissipation  functions 
to  achieve  numerical  stability.  The  resulting  set  of  equations  is  solved  using  an  implicit, 
factored,  time-stepping  algorithm.  The  solution  takes  place  on  a  computational  grid  that 
is  generated  around  the  configuration  in  zones;  zonal  boundaries  are  transparent  to  the 
flowfield. 

The  ram  accelerator  projectile  geometry  shown  in  Figure  2  represents  a  simplification 
of  the  actually  configuration  used  for  test  firings.  The  actual  projectile  includes  a  set  of 
four  bore-riding  fins  that  extend  from  the  point  of  maximum  projectile  diameter  to  the 
projectile  base  and  span  the  area  between  the  projectile  and  the  launch  tube  (see  Figure  1). 
Exclusion  of  these  fins  permits  a  2D/axisymmetric  calculation  and  results  in  a.  significant 
computer  time  savings  over  the  full  3D  numerical  simulation.  Such  3D  simulations  are 
possible  with  the  present  code  (e.g.  see  Reference  9);  application  to  the  ram  accelerator 
will  be  addressed  in  future  reports. 


1.  Equations  of  Motion. 


The  RANS  equations  for  2D/axisymmetric  flow  are  written  in  the  following  conser¬ 
vation  form.  The  dependent  variables  u,  v  and  e  are  mass- averaged. 


dW  dF  dG  (G 

dt  ***  dx  dy  \  y 


(1) 


16  Ota,  D.K.,  Chakravarthy  S.R.,  and  Darling,  J.C.,  ‘Mn  Equilibrium  Air  Navier-Stokes  Code  for  Hypersonic  Hows,"  AIAA- 
88-0419,  Proceedings  of  the  26th  AIAA  Aerospace  Sciences  Meeting ,  Renp,  NV,  January  11-14,  1988. 

1 7  Palaniswamy,  S.,  and  Chakravarthy,  S.R.,  “Finite  Rate  Chemistry  for  USA  Series  Codes:  Formulation  and  Applications 
AlA  A-89-0200,  Proceedings  of  the  27th  AIAA  Aerospace  Sciences  Meeting,  Reno.  NV,  January  9-12,  1989. 

18  Baldwin,  B.S.  and  Lomax,  II.,  “ Thin  Layer  Approximation  and  Algebraic  Model  for  Separated  Turbulent  Flows,”  AIAA- 
78-257,  Proceedings  of  the  16th  AIAA  Aerospace  Sciences  Meeting,  Huntsville,  AL,  January  16-18,  1978. 

19 Goldberg  U.C. ,  “ Separated  Flow  Treatment  with  a  New  Turbulence  Model ,,f  AIAA  Journal.  Vol.  24,  No.  10,  October 
1986,  pp.  1711-1713. 
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where  a  =  1  for  axisymmetric  flow  and  0  for  two-dimensional  flow. 
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In  Equations  2-7,  the  laminar  and  eddy  viscosities,  p  and  pt,  are  implicitly  divided  by 
the  reference  Reynolds  number.  The  equations  used  for  the  Euler  (inviscid)  calculations 
are  obtained  from  Equations  1-9  by  setting  both  laminar  and  eddy  viscosities  to  zero.  The 
flow  medium  is  assumed  to  be  a  perfect  gas  satisfying  the  equation  of  state, 


p  =  p3?T 


(10) 


The  following  power  law  was  used  to  relate  molecular  viscosity  to  temperature:  20 


P_ 

Po 


(ID 


where  p0  =  0.1716  mP,  T0  =  491.6  R,  and  n  =  0.64874.  The  laminar  and  turbulent  Prandtl 
numbers,  Pr  and  Pr(,  were  assumed  constant  with  values  of  0.72  and  0.9  respectively.  The 
ratio  of  specific  heats,  7,  was  also  assumed  constant.  The  specific  heat  capacities  at 
constant  volume  and  pressure,  cv  and  cp,  are  related  as  7  =  cp/cv  and  3?  =  (7  -  l)cp/7- 

20  Wa zor  G.,  Ben-Dor  G.,  and  Igra  O.,  “A  Simple  uni  Accurate  Expression  for  the  Viscosity  of  Nonpolar  Diatomic  Gases 
up  to  10,000  K."  A I A  A  Journal.  Vol.  S3,  No.  4,  April  1985,  pp.  636-638. 
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Assuming  a  time-invariant  grid  and  using  the  transformation  of  coordinates  implied 
by  r  =  t,£  =  ((x,y),rj  =  rj(x,y),  Equation  1  can  be  recast  into  conservation  form  where  £ 
and  r/  are  the  new  independent  variables  and  x$,  xv ,  and  yv  are  the  four  transformation 
coefficients  obtained  numerically  from  the  mapping  procedure.8 

fr  +  +  <-<  +  +  f  -  f)  =  0  (12) 

The  “Area”  in  Equation  12  denotes  the  area  of  the  finite  volume  cell  under  consideration 
at  the  time  of  discretization  of  the  equations.  The  transformed  time  variable  is  represented 
by  r. 


2.  Turbulence  Modeling. 

The  shock/boundary-layer  interference  flowfield  between  projectile  and  launch  tube 
as  well  as  the  projectile  wake  can  include  regions  of  recirculating  flow.  Modeling  of  these 
regions  can  be  critical  to  the  overall  flowfield  solution  quality.  However,  most  existing 
turbulence  models  either  do  not  treat  such  regions  or  do  so  in  a  semi-empirical  fashion 
that  is  frequently  inadequate.  To  improve  the  predictive  capability  of  separated  flows 
using  RANS  codes  a  new  turbulence  model  has  been  recently  developed.19  The  new  model 
is  based  on  experimental  observations  of  detached  flows.  The  model  prescribes  turbulence 
kinetic  energy  ( k )  and  dissipation  (e)  analytically  within  backflows.  A  Gaussian  variation 
of  k  normal  to  wall  surfaces  is  assumed.  The  length  scale  of  turbulence  is  proportional  to 
the  local  distance  from  the  wall  to  the  edge  of  the  viscous  sublayer  which  is  located  outside 
the  backflow  region.  The  latter  feature  is  a  basic  assumption  of  the  model.  The  stress 
scale  is  the  local  maximum  Reynolds  stress  which  typically  occurs  around  the  middle  of 
the  boundary  layer  well  outside  the  separation  bubble.  This  scale  must  be  supplied  by  a 
turbulence  model  that  is  used  beyond  backflow  regions. 

The  main  equations  of  the  backflow  model  are  given  in  Reference  19.  A  formula  for 
the  eddy  viscosity  (/z*)  distribution  within  backflows  is  derived  and  used  for  the  RANS 
equations  when  the  calculations  are  done  inside  separation  bubbles.  Outside  of  them, 
another  turbulence  model  (for  example  Baldwin-Lomax18)  supplies  the  values  of  eddy 
viscosity.  While  the  Baldwin-Lomax  turbulence  model  is  used  to  detect  flow  separation  and 
to  initiate  application  of  the  backflow  model,  the  latter  model  can  relocate  the  separation 
point.  Comparisons  of  the  Baldwin-Lomax  turbulence  model  and  Goldberg's  backflow 
model  are  given  in  References  10  and  21.  For  further  details  of  how  the  model  treats 
the  influence  of  large  eddies  residing  outside  detached  regions,  the  history  effect  of  these 
eddies  downstream  of  reattachment,  and  the  mutual  influence  of  multiple  walls  on  the  eddy 
viscosity,  see  References  19  and  22. 

31  Goldberg,  U.C.,  “ Separated  Flows  Calculations  Wi'lA  A  New  Turbulence  Model,"  presented  at  the  IACM  First  World 
Congress  on  Computational  Mechanics,  TX,  Sept,  1986. 

33  Goldberg,  U.C.,  “ Prediction  0}  Separated  Flows  With  A  New  Turbulence  Model.”  AIAA  Journal.  Vol.  26,  No.  i.  April 
1988,  pp.  405-408. 


6 


3.  Computational  Algorithm. 

The  spatial  discretization  technique  for  the  equations  of  motion  (Equation  1)  must  be 
reliable  and  robust  if  it  is  to  successfully  capture  the  complex  physics  of  projectile/launch 
tube  interacting  flowfields.  The  TVD  formulation  for  the  convective  terms  (the  hyperbolic 
part  of  the  time-dependent  Navier-Stokes  equations)  along  with  a  special  treatment  of  the 
dissipative  terms  provides  an  appropriate  simulation.  In  recent  years,  TVD  formulations 
have  been  constructed  for  shock-capturing  finite-difference  methods.8,9  Near  large  gradi¬ 
ents  in  the  solution  (extrema)  TVD  schemes  automatically  reduce  to  first-order  accurate 
discretizations  locally  while  in  general  they  can  be  constructed  to  be  of  higher-order  accu¬ 
racy.  This  local  effect  restricts  the  maximum  global  accuracy  possible  for  TVD  schemes 
to  third  order  for  steady-state  solutions. 

TVD  methods  manifest  many  properties  desirable  in  numerical  solution  procedures. 
By  design  they  avoid  numerical  oscillations  and  “expansion  shocks”  while  at  the  same  time 
being  of  higher-order  accuracy;  “expansion  shocks”  are  shock  waves  which  do  not  satisfy 
the  entropy  inequality.  TVD  formulations  are  also  based  on  the  principle  of  discrete  or 
numerical  conservation  which  is  the  numerical  analog  of  physical  conservation  of  mass, 
momentum,  and  energy.  Thus,  TVD  schemes  can  “capture”  flowfield  discontinuities  (e.g. 
shock  waves)  with  high  resolution.  At  a  fundamental  level  they  are  based  on  upwind 
schemes;  therefore,  they  closely  simulate  the  signal  propagation  properties  of  hyperbolic 
equations.  Schemes  based  on  the  TVD  formulation  are  completely  defined.  In  contrast, 
schemes  based  solely  on  central  differences  involve  global  dissipation  terms  for  stability  and 
have  one  or  more  coefficients  that  must  be  judiciously  chosen  to  achieve  desirable  results. 

Any  conventional  time  discretization  method  suitable  for  the  Navier-Stokes  equations 
can  be  used  together  with  this  space  discretization  methodology;  for  example,  approxi¬ 
mate  factorization  and  relaxation  techniques.  Proper  treatment  of  the  dissipative  terms 
of  the  Navier-Stokes  equations  is  also  important  in  the  construction  of  reliable  numerical 
methods.  Unidirectional  second  derivative  terms  are  treated  by  using  central  difference 
approximations.  Cross  derivatives  are  represented  by  finite-differences  the  nature  of  which 
depends  upon  the  sign  of  the  coefficient  of  such  terms.  This  treatment  augments  diagonal 
dominance  of  the  resulting  set  of  discretized  equations  without  detracting  from  the  accu¬ 
racy  and  while  adding  to  the  reliability  of  the  numerical  procedure.  Further  details  can  be 
found  in  References  8  and  9. 


4.  Computational  Grid. 

Computing  in-bore  projectile  flowfields  is  complicated  by  the  multi-wall  geometry.  As 
illustrated  in  Figures  1  and  2  the  ram  accelerator  projectile  consists  of  several  sharp  corners 
that  would  severely  hamper  conventional  grid  generation  schemes  that  require  one  set  of 
grid  lines  to  be  tangent  to  surfaces  and  another  set  to  be  normal  to  them.  This  geometry  is 
more  easily  gridded  by  the  zonal  approach.  The  internal  geometry  of  the  ram  accelerator 
launch  system  in  broken  up  into  two  zones  of  simple  geometric  shape  (Figure  4).  In  each 
zone  any  approach  to  generating  the  grid  can  be  used  such  as  algebraic  methods,  differential 
equations  methods,  etc.  In  this  case  an  algebraic  grid  is  used  with  clustering  near  surfaces 
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and  other  regions  in  the  flowfield  where  high  gradients  are  expected  (Figure  5).  In  the 
zonal  approach,  the  computational  method  and  computer  program  are  constructed  in  such 
a  manner  that  each  zone  may  be  considered  as  an  independent  module  interacting  with  each 
other  before  or  after  the  information  corresponding  to  each  zone  is  updated  one  cycle.  In 
addition,  the  zonal  boundaries  are  transparent  to  flowfield  phenomena  (e.g.  shock  waves). 
The  actual  grid  used  for  these  computations  (~  30,000  nodes)  consisted  of  the  following 
dimensions  for  zones  1  and  2:  457x51,  201x21.  The  computational  domain  starts  a  small 
distance  forward  of  the  body  and  extends  1.5  projectile  body  lengths  downstream.  The 
computational  results  were  found  to  be  essentially  independent  of  further  grid  refinement 
when  all  other  factors  were  the  same.  A  dense  grid  is  preferred  for  resolving  flow  details 
within  boundary  layers,  to  prevent  shock  smearing,  and  to  resolve  shock/boundary-layer 
interactions.  Grid  zones  have  been  generated  in  such  a  manner  that  the  vertical  families 
of  mesh  lines  are  continuous  between  zones  1  and  2  aft  of  the  projectile.  While  this  type 
of  grid  is  natural  for  this  problem  the  computational  methodology  permits  a  wider  class 
of  patched  zonal  grids  where  neither  horizontal  nor  vertical  grid  lines  need  be  continuous 
across  zones. 


5.  Results. 

Numerical  simulations  for  non-reacting  flow  were  performed  using  the  four  geometries 
listed  with  Figure  2.  The  diffuser  area  ratio  is  defined  as, 


Diffuser  Area  Ratio  = 


_ (tube  diameter)2 _ 

(tube  diameter)2  — (max.  projectile  diameter)2 


Small  ratios  indicate  large  clearance  between  projectile  and  tube  walls.  Calculations  were 
performed  on  the  BRL  CRAY-2  supercomputer  using  1.0  to  1.3  CPU  hours  per  steady- 
state  solution.  Upstream  boundary  conditions  are  specified  (slug  flow)  while  supersonic 
outflow  conditions  are  extrapolated.  Numerical  simulations  were  done  in  the  projectile- 
fixed  reference  frame  thus,  a  no-slip  condition  is  employed  on  the  projectile  body  while  the 
tube  wall  translates  at  the  freestream  velocity.  The  plotted  results  have  been  transformed 
to  the  inertial  reference  frame. 

Figure  6  shows  the  computed  pressure  contours  for  the  38mm  configuration  1,  ve¬ 
locity  of  1102  m/s  (M  —  3.022),  and  tube  fill  pressure  of  25  atm  (Re  =  2.3xl08  based 
on  body  length).  An  oblique  shock  is  attached  to  the  nose  of  the  projectile.  This  shock 
reflects  from  the  tube  wall  and  again  from  the  turbulent  boundary  layer  on  the  body.  Mul¬ 
tiple  shock/boundary-layer  interactions  occur  downstream.  Contours  in  the  wake  show 
a  weakening  of  shock  reflections.  The  horizontal  line  extending  from  the  projectile  base 
downstream  is  the  zonal  grid  boundary  and  is  transparent  to  the  flowfield.  Figures  7  and  S 
show  the  computed  pressure  contours  for  the  38mm  configuration  2  for  a  velocity  of  1102 
m/s  ( M  =  3.022)  tube  pressure  25atm,  and  1430  m/s  ( M  =  3.95)  tube  pressure  20atm, 
respectively. 

Figure  9  shows  the  computed  tube  wall  pressure  for  the  38mm  configuration  2,  M 
=  3.022  and  tube  fill  pressure  25  atm.  Computations  for  inviscid,  laminar,  and  turbulent 
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flow  are  compared  to  pressures  measured  at  the  University  of  Washington  (unpublished). 
Overall  the  comparison  is  quite  good  even  in  the  wake  region  (x  >  .17m).  The  pressure 
rise  between  .05  <  x  <  .075m,  caused  by  the  initial  reflection  of  the  body  nose  shock  on 
the  tube  wall,  is  simulated  but  at  a  reduced  magnitude  (perhaps  due  to  the  steady  flow 
assumption).  Two  pressure  peaks  between  .075  <  x  <  .lm  are  computed  with  the  second 
slightly  overpredicting  the  measured  peak  at  x  ~  .lm.  These  peaks  arise  from  secondary 
reflection  pf  the  body  nose  shock  on  the  tube  wall.  The  simulations  produce  noticeable 
but  insignificant  differences  in  wall  pressure  for  inviscid,  laminar,  and  turbulent  modeling; 
however,  projectile  surface  and  tube  centerline  (wake)  pressures  are  affected  to  a  greater 
degree  by  boundary  layer  modeling  as  shown  in  Figures  10  and  11.  Pressure  oscillations 
on  the  body  surface  (Fig.  10,  0  <  x  <  .07m)  are  caused  by  pressure  signals  that  travel  up¬ 
stream  in  the  near-wall  flowfield  and  may  be  related  to  the  use  of  a  body-fixed  coordinate 
system.  Figure  12  shows  tube  wall  pressure  distributions  for  both  38mm  configurations. 
Although  the  general  trends  are  the  same,  the  pressure  levels  for  configuration  2  are  sig¬ 
nificantly  larger.  Configuration  2  has  a  18%  larger  diffuser  area  ratio  (i.e.  smaller  diffuser 
area  resulting  in  higher  pressures)  and  boundary  layer  thicknesses  are  different  as  well. 
Figure  13  shows  tube  wall  pressure  for  M  =  3.022  and  3.95.  At  the  higher  velocity  the 
body  nose  shock  intersects  the  tube  wall  farther  downstream  and  causes  shock  reflections 
that  result  in  lower  pressures  along  the  tube  wall.  Low  Mach  number  (M  =  3.022)  simula¬ 
tions  for  configuration  2  produce  a  double  pressure  pulse  at  x  ~  .lm;  this  is  caused  by  the 
intersection  of  a  reflected  shock  from  the  tube  wall  near  the  end  of  the  projectile  forebody 
(see  Figure  7).  In  the  wake  region  (x  >  .17m)  tube  wall  pressures  are  nearly  invariant  for 
these  Mach  numbers. 

Figure  14  shows  tube  wall  pressure  comparisons  for  the  120mm  configuration  2  using 
inviscid  and  turbulent  flow  and  the  same  flow  conditions  as  previous  figures  (Re  =  9.7xl08 
based  on  body  length).  Inviscid  flow  modeling  yields  higher  pressures  downstream  of 
initial  and  secondary  shock  reflections  (x  >  ,3m)  but  lower  pressures  in  the  wake  region 
(x  >  .54m).  Figure  15  shows  the  comparison  of  tube  wall  pressures  for  both  configurations 
of  the  120mm  geometry.  Configuration  2  shows  slightly  higher  pressures  than  configuration 
1  since  the  diffuser  area  is  15%  larger.  A  similar  comparison  for  the  38mm  geometry 
(Figure  12)  showed  large  differences  in  tube  wall  pressures.  The  projectile  forebody  and 
afterbody  angles  are  similar  for  120mm  configurations  but  noticeably  different  for  38mm 
configurations.  Figure  16  shows  120mm  configuration  2  tube  wall  pressure  for  M  =  3.022 
and  3.95.  At  the  higher  velocity  the  body  nose  shock  intersects  the  tube  wall  farther 
downstream  but  causes  shock  reflections  that  result  in  similar  pressures  along  the  tube 
wall.  Figure  17  shows  tube  wall  pressures  for  38mm  configuration  2  (compared  to  data 
in  Fig.  9)  and  120mm  configurations  1  and  2  (x  non-dimensionalized  by  the  total  body 
length,  L ).  Due  to  different  scales  for  the  38mm  and  120mm  geometries,  Reynolds  number 
effects  are  significant.  Smaller  differences  in  tube  wall  pressure  for  38mm  and  120inm 
geometries  result  from  higher  Mach  number  (M  =  3.95)  simulations  as  illustrated  in  Figure 
18;  however,  large  differences  are  apparent  in  the  wake  flow  ( x/L  >  1)  for  both  Mach 
number  cases. 
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IV.  REACTING  FLOW 


The  BRL-LFD  RAMCOMB  (RAMjet  COMBustion)  code  has  been  used  for  the  nu¬ 
merical  simulation  of  solid  fuel  combustion  in  a  tubular  solid-fuel  ramjet  (SFRJ)  projectile.11 
Solid  fuel  regression  rate  and  projectile  thrust  predictions  compare  favorably  with  in-flight 
and  ground  test  data.  For  the  SFRJ  application  the  RAMCOMB  code  simulated  a  mass- 
controlled  (stoichiometric)  reaction  of  non-premixed  solid  fuel  and  oxygen  using  classical 
diffusion  flame  assumptions.  A  reaction  rate  expression  was  not  required  since  the  reac¬ 
tion  was  diffusion-controlled  and  the  reactants  were  fully  consumed  and  thus,  could  not 
coexist  at  any  point  in  the  flowfield.  For  application  to  the  ram  accelerator  system  the 
RAMCOMB  code  was  modified  for  premixed,  finite-rate,  gaseous  fuel  combustion.  Reac¬ 
tion  rates  are  formulated  in  terms  of  temperature  and  species  mass  fraction.  Reactants 
and  products  may  coexist  and  need  not  be  fully  consumed  in  the  reaction.  RAMCOMB 
simulates  subsonic  flow  and  thus,  reacting  flow  simulations  are  constrained  to  the  flowfield 
downstream  of  the  normal  shock  (Figure  3). 


1.  Equations  of  Motion. 


The  ram  accelerator  projectile  configuration  can  support  a  region  of  low  speed  recir¬ 
culating  flow  in  the  subsonic  flowfield  downstream  of  the  normal  shock  (Figure  3).  Thus 
the  conservation  equations  are  elliptic  in  form.  Since  the  geometry  is  axisymmetric  the 
equations  can  be  written  in  cylindrical  coordinates.  The  velocity  components  in  this  sys¬ 
tem  are  u,u,  and  w  for  the  radial  (r),  azimuthal  ( 9 ),  and  axial  (z)  directions,  respectively. 
Axisymmetric  flow  is  assumed  thus,  all  ^-derivatives  are  ignored;  however,  the  azimuthal 
velocity  component  and  the  azimuthal  momentum  equation  are  retained.  Since  steady 
flows  are  considered,  time  derivatives  ( d/dt )  are  ignored.  The  conservation  equations  for 
mass  (global)  and  momentum  are  given  by,23 


V-pV-  +  d(pw)  =0 

r  or  oz 

(13) 

V  •  [puV  —  Tr]  —  ~[pv2  -  Tgg]  +  ^  =  0  (radial) 

r  Or 

(14) 

V  •  [pwV  —  tz]  +  ~  =  0  (axial) 

Oz 

do) 

-V  •  [r{ovV  —  f^)]  =0  (azimuthal) 

1" 

(1C) 

Energy  conservation  for  a  compressible  flow  is  expressed  by  the  First  Law  of  Thermo¬ 
dynamics.  The  steady  form  of  the  First  Law  states  that  the  net  rate  of  stagnation  enthalpy 
( h )  inflow  for  a  control  volume  is  equal  to  the  sum  of  the  shear  work  done  by  the  contents 
of  the  control  volume  on  the  surroundings  (r)  and  the  heat  transfer  to  the  surroundings.23 


V 


pVh  +  Jh  +  52  hjJj  +  Jk  —  (ufr  +  VTq  +  WTZ) 


23 Schlichting,  //.,  Boundary  Layer  Theory,  7th  ed.,  translated  by  J.  Kestm,  McGraw-Hill,  New  York,  1979. 
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where  mass,  heat,  and  turbulence  kinetic  energy  flux  are  defined  as, 

Jj  =  (/Xeff/R  e)Vm;  =  TjVmj 
Jh  =  (^eff/Pr)cP;  VT  =  ThcP]VT 
Jk  =  (pca/Pv)Vk  =  TkVk 

The  mass  fraction  and  molar  specific  enthalpy  for  species  j  are  mj  and  hj,  respectively. 
Radiation  flux  is  neglected  in  Equation  17,  but  will  be  explored  in  future  reports. 

In  Equations  14-17  the  shear  stress  (r)  includes  the  Reynolds  stress  with  an  effective 
fluid  viscosity  expressed  as  the  sum  of  the  molecular  and  turbulent  viscosities,  /ieff  = 
/x  +  fJLt .  The  subscript  “eff”  serves  as  a  reminder  that  the  turbulent  fluid  properties  are 
being  treated  as  kinematic  properties  with  a  set  of  augmented  transport  coefficients.  The 
calorically  perfect  gas  assumption  is  made  since  the  temperature  dependence  of  cP]  for  the 
reactants  and  products  is  not  well  determined.  The  Prandtl  number  (Pr)  is  assumed  to  be 
nearly  unity  (.9)  which  is  considered  adequate  for  gaseous  flows  even  with  combustion21. 


For  flows  that  are  subsonic  and  dominated  by  recirculation  the  stream  function- 
vorticity  form  of  the  governing  equations  has  been  widely  utilized.23  This  form  facilitates 
the  use  of  numerically  efficient  Gauss-Seidel  relaxation  algorithms.  Stream  function,  v. 
and  vorticity,  u  are  defined  using,23 

dx b 

—  =  rpw  (IS) 

or 


—  =  —rpu 


_  du 

UJ  =  V  X  V  =  - - 

dz 


1  dxp 
rp  dz 


d_ 

dr 


1  dtp 

rp  dr 


The  governing  equations  are  derived  in  Reference  11  and  can  be  expressed  in  the  form 
of  a  general  variable  <f>.  This  variable  can  represent  stream-function,  vorticity,  azimuthal 
velocity,  stagnation  enthalpy,  or  species  mass  fraction. 


d  ( Ad^\  d  (a^W  5  [  d  1  d  \k  d  ( 
a*  b*rTzM)  hra-r(c' 


+  rd j,  =  0 


II 

a0  =  0, 

b*  =  1  /(pr2), 

C<t>  ^  1 

dj,  =  —uj/r 

s-T 

II 

a4>  =  r2, 

h  =  r2, 

C<t>  /^efF  1 

dip  Qui 

<A  =  rv, 

® <t>  1  > 

^ 4>  r  , 

c <t>  =  1/r2, 

dj  =  0 

»-cf 

II 

—  I? 

b<)>  (Fh)eff  » 

c<t> 

d*  =  Qh  —  Si 

5 

II 

I7 

(  Tj  )eff  i 

c<i> 

1 

II 

24  Bradshaw,  P.,  Cebeci,  T.,  and  Whitelaw,  J.H. ,  Engineering  Calculation  Methods  for  Turbulent  Flows,  Academic  Press, 
New  York,  1981. 
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is  defined  in  Reference  11. 


o  =  A<eff/r 

i?j  =  reaction  rate  per  unit  volume  for  species  j. 


For  N  species  only  N  —  1  specie  equations  (<j>  =  rrij )  must  be  solved,  since  the  sum  of  the 
mass  fractions  must  equal  unity.  In  effect  the  continuity  equation  (Eq.  13)  is  the  jVth 
specie  equation  since  the  summation  of  all  specie  equations  yields  the  continuity  equation. 
The  mixture  equation  of  state  follows  from  Dalton’s  Law, 


V  =  P%TY. 

i 


(22) 


where  &  =  dZEjMj,  Mj  is  the  molecular  weight  of  species  j,  and  3?  is  the  specific  gas 
constant.  Mixture  temperature  (T)  is  obtained  from  the  conservation  of  energy  expressed 
in  terms  of  the  stagnation  enthalpy, 


h  =  T^cPjm7  + 


V2  1 
2  +  l  Sc 


£! 

PrJ  2  + 


1_ 

LSc 


l 

Pr 


,m  j 


(23) 


with  V  =  (u2  +  v2  +  w2y'2,  and  V  is  the  magnitude  of  the  turbulent  (fluctuating)  velocity. 
The  Schmidt  number  (Sc)  is  assumed  to  be  unity.  Mixture  viscosity  (p)  is  defined  using 
Sutherland’s  expression25  for  T  <  3400R  and  using  Equation  11  for  higher  temperatures. 


A* 


2.270xl0“8 


j'l.S 

T+  198.6 


(24) 


2.  Turbulence  Modeling. 

A  two-equation  turbulence  model  suggested  by  Kim  and  Chung26  for  reacting  flows  has 
been  added  to  the  code.  This  model  describes  the  turbulence  viscosity  (/zt)  as  a  function 
of  turbulence  kinetic  energy  ( k )  and  dissipation  rate  (e)  variables.  A  set  of  partial  differ¬ 
ential  equations  is  written  for  k  and  t  and  solved  along  with  the  Navier-Stokes  equations 
(Equation  21). 


dk 


dk  1 


pWTz+f,UTr~r 


d_ 

dz 


=  G  —  pe 


(25) 


25/tme«  Research  Staff,  “ Equations ,  Tables,  and  Charts  for  Compressible  Flow,*  SAC  A  Report  1135,  195$. 

26  Aim,  Y.M.,  and  Chung,  T.J.,  “Finite- Element  Analysts  of  Turbulent  Diffusion  Flames.”  A!  A  A  Journal,  IW.  27,  .Vo.  3. 
March  1989,  pp.  330-339. 
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where,  p,k  =  p  +  fit/ h,  =  fi  +  fit/ KAk  =  l,At  =  1.3,  C\  =  1.44,  C2  =  1.92,  C3  =  .09. 

Results  using  this  turbulence  model  will  included  in  future  reports. 


3.  Combustion  Modeling. 

Systems  of  chemical  reactions  are  usually  complex.  A  hydrocarbon  fuel,  for  example, 
may  contain  hundreds  of  distinct  chemical  species.  During  combustion  numerous  short¬ 
lived  intermediate  species  axe  formed  in  addition  to  the  final  products  that  include  H20, 
CO,  and  C02.  Inclusion  of  intermediate  species  and  reaction  steps  within  the  framework 
of  a  calculation  could  be  defeated  by  an  ignorance  of  fhe  thermodynamic,  transport,  and 
chemical-kinetic  properties  of  these  species.  It  is  possible  to  make  useful  predictions  of 
combustion  processes  by  confining  attention  to  a  postulated  global  reaction  between  fuel 
and  oxidizer  that  neglects  intermediate  steps. 

Fuel  -f  ^Oxidizer  +  Diluents  =>  (1  4-  (^Products  4-  Diluents  (29) 

where  S  is  the  stoichiometric  oxidizer/fuel  ratio.  Two  limitations  of  a  global  reaction 
model  should  be  noted.27  During  the  ignition  delay  period,  when  the  complete  reaction 
mechanism  will  predict  the  gradual  buildup  of  free  radicals  with  little  or  no  perceptible 
temperature  change,  a  global  model  will  indicate  an  immediate  exothermic  (or  endother¬ 
mic)  reaction  with  resulting  temperature  increase  (or  decrease).  Secondly,  the  adiabatic 
flame  temperature  predicted  by  the  global  model  will  be  higher  than  for  the  complete 
multi-step  reaction  mechanism.  This  results  from  the  fact  that  a  global  model  does  not 
include  all  of  the  product  species  actually  present  in  the  reaction.  If  very  low  pressures 
are  avoided,  the  discrepancy  in  flame  temperature  will  not  be  significant. 

Many  of  the  assumptions  commonly  made  in  combustion  models28-30  are  also  incor¬ 
porated  into  the  present  model;  in  particular,  assumptions  made  in  the  Shvab-Zcldovich 
formulation.30  The  reacting  fluid  is  treated  as  a  continuum  with  equal  mass  diffusivity 
properties  for  all  species.  The  specific  heat  at  constant  pressure  (cPj)  of  each  species  is 
assumed  constant  in  the  gas  phase  (see  Equation  23).  The  contribution  of  turbulence  is 
included  as  turbulence  kinetic  energy  (Equation  21,  <t>  =  k)  and  the  fluctuating  velocity 
component  (F)  in  Equation  23.  A  one-step,  irreversible,  temperature  and  mass  controlled, 
finite-rate  reaction  is  considered  for  which  the  reactant  species  need  not  be  fully  consumed. 
The  reactants  are  gaseous  fuel  (CH4)  and  oxidizer  (02)  premixt~i  to  stoichiometric  pro¬ 
portions.  The  products  of  this  reaction  (e.g.  C02  and  H20)  are  lumped  into  one  specie 

27  Rogers,  R.C.,  and  Chmitz,  W.,  “Uaing  a  Global  Hydrogen- Air  Combustion  Model  in  Turbulent  Reacting  Flow  Calcula¬ 
tions,'*  AJ_A_A_JournaL,  Vol.  21,  So.  4,  April  1988,  pp.  586-592. 
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with  averaged  thermodynamic  properties  ( M.  determined  by  a  weighted  average).  Diluent 
species  are  included  in  Equations  22  and  23. 

Hydrocarbon  reactions  are  commonly  used  for  ram  accelerator  testing  at  the  Univer 
sity  of  Washington,1  for  example, 

-CH4  +  202  +  6N2  =*  C02  +  2H20  +  6N2  +  7C2  +  3H2  (30) 

2  4 

This  reaction  is  one  of  a  general  class  of  hydrocarbon  reactions  investigated  by  Westbrook 
and  Dryer,31 

CaHbOc  4*  ~ (a  ■+■  2  +  c)02  =>  aCO  +  — H20  (31 ) 

where  elemental  and  diluent  species  N2,  C2  and  H2  are  ignored  thus, 

CH4  +  202  =»  C02  +  2H20  (32) 

CH4  is  considered  the  fuel,  02  the  oxidizer,  and  products  are  C02  and  H20.  Thermody¬ 
namic  data  is  given  in  Table  1. 


Table  1.  Thermodynamic  Data 


Cp(  mnl*  K  ) 

M 

ch4 

8.536 

-17895 

16.043 

02 

7.017 

0 

32.000 

co2 

8.874 

-94054 

44.011 

h2o 

8.025 

-57798 

18.016 

Fuel 

8.536 

-17895 

16.043 

Oxidizer 

7.017 

0 

32.000 

Product 

8.450 

-75926 

26.681 

Elemental  species  (02,  N2,  C2,  H2)  do  not  contribute  to  the  heat  of  reaction  since 
A H°j  =  0  but  are  included  in  the  equations  of  state  (Equations  22  and  23).  The  reaction 
is  exothermic  with  a  heat  of  reaction  given  by,29 


AHr 


(x>^n)  - 

V  3  /  prod  V  i  /  i 


(33) 


where  n3  is  the  stoichiometric  coefficient  and  A HJ  is  the  heat  of  formation  for  each  species 
j  in  Equation  30.  For  this  reaction  AHr  =  -690224.8  J.  For  a  methane/oxygen  mixture 

20  Khalil,  E.E.,  Spalding,  D.B. ,  and  Whitclaw,  J.H.,  “The  Calculation  of  Local  Flow  Properties  in  Two-Dimensional  Fur- 
naces,  **  International  Journal  of  Heat  and  Mass  Transfer ,  1975,  Voi  IS,  pp.  775-791. 

*9  Kuo,  K.K.,  Principles  of  Combustion,  John  Wiley  and  Sons,  Inc.,  New  York,  1986. 

30  Williams ,  F.A.,  Combustion  Theory,  Addison- Wesley  Publishing  Company,  Inc.,  1965. 

31  Westbrook,  C.K.,  and  Dryer,  F.L.,  u Chemical  Kinetic  Modeling  of  Hydrocarbon  Combustion ,*  Progress  in  Energy 
Combustion  Science.  Vol.  10,  pp.  1-57. 
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at  300°K,  the  heat  of  reaction  is  nearly  invariant  with  pressure  above  10  atm.32  The 
combustion  flame  temperature  is  given  by,29 


AH'  +  ~  hT,)j  I  =0  (34) 

'  j  /  prod 

where  Tj  is  the  adiabatic  frozen-flame  temperature  and  TJ  is  the  initial  or  base-state 
temperature  (e.g.  298K).  For  the  particular  reaction  (Equation  30)  Tj  can  be  iteratively 
computed  to  be  1643. 06K  (2957.50R). 


Conservation  of  chemical  species  can  be  written  for  each  specie  in  Equation  32.  For 
example,  the  conservation  of  mass  for  product  CO2  is  given  by  (see  Equation  21 ,  4>  =  mco2  )< 


d 

dtp 

d 

dtp 

d 

rn  ,  dmco2 

dz 

mc°j 

dr 

mc0la7 

dz 

(ico  2)effr 

(rCo2) 


eff  ' 


dmc  o2 
dr 


—  rRcOi  =  0 


(35) 

where  T  is  a  diffusion  coefficient  (/xeff/ Re)  and  i?co2  is  the  reaction  rate  per  unit  volume. 
The  reaction  rate  is  defined  using  the  Law  of  Mass  Action,* 


N 

R—C  n  m] ’ 

1 


(36) 


where  C  is  the  specific  reaction  rate  constant  and  iV  =  4  for  Equation  32.  Using  an 
Arrhenius  expression  for  C  the  reaction  rate  is  approximated,31 


R  =  ATaex  p 


m 


a 

CH< 


(?7) 


where  ATa  is  the  collision  frequency,  the  exponential  term  is  the  Boltzmann  factor,  and 
Ea  is  the  activation  energy.  For  the  general  hydrocarbon  reaction,  Reference  31  specifics 
A  =  8.3xl05  moles/cm3-sec,  a  =  0,  Ea  =  125520  J/mole,  a  =  -.3,  b  =  1.3,  and  -ft  =  8.314 
J/K-mole.  The  magnitude  of  the  collision  frequency  determines  the  speed  at  which  the 
reaction  progresses. 


4.  Computational  Algorithm. 

Equation  21  can  be  reduced  to  a  successive-substitution  formula  for  flow  variable 
4>'  at  each  node  on  the  computational  grid.  Central  finite  differences  are  used  for  the 
diffusive  and  source  terms  and  upwind  differences  for  the  convective  terms.  Using  upwind 
differencing  in  the  specie  conservation  equations  reduces  the  occurrence  of  negative  specie 
mass  fractions  in  mixing  layers.  The  resulting  system  of  equations  for  the  entire  grid  is 
solved  using  a  Gauss-Seidel  relaxation  scheme.11  Each  iteration  cycle  is  made  up  of  A 
sub-cycles,  where  A  is  the  number  of  equations  being  considered  (A  must  be  at  least  2  since 
the  equations  for  4>  =  u/r  and  $  =  xj>  are  the  minimum  required  to  define  the  flow).  In 

33  Kivtty,  Y - ,  Grimbtrg,  <Vf.,  and  Halevy,  D.,  “Performance  of  the  Ram  Accelerator  in  the  Oblique  Detonation  Mode," 
Proceeding »  of  the  flit  Meeting  of  the  Aeroballistic  Range  Association,  San  Diego,  CA,  Oct.  ££-25,  1990. 

*  “The  rate  of  disappearance  of  a  chemical  species  is  proportional  to  the  products  of  the  concentrations  of  the  reacting 
chemical  species,  each  concentration  being  raised  to  a  power  equal  to  the  corresponding  stoichiometric  coefficient"39 
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each  sub-cycle  grid  points  are  scanned  row  by  row  and  a  single  variable  is  updated.  The 
variables  ui/r  and  rp  are  updated  in  order  followed  by  all  other  variables.  When  all  sub¬ 
cycles  are  completed  a  new  iteration  cycle  is  started  in  which  the  values  of  the  variables 
from  the  latest  iteration  are  immediately  used.  This  is  consistent  with  the  Gauss-Seidel 
methodology.33  Convergence  is  satisfied  when  the  greatest  relative  change  in  any  flow 
variable,  <f>,  over  all  grid  points  is  smaller  than  a  prescribed  tolerance. 


5.  Results. 

Reacting  flow  simulations  were  performed  for  two  geometries,  the  38mm  configuration 
2  and  the  120mm  configuration  2  (see  Figure  2).  The  reaction  considered  is  given  by 
Equation  30  and  the  flow  conditions  were  1430  m/s  velocity,  20  atm  pressure  and  M 
=  3.95.  For  numerical  simulations  the  mixture  is  assumed  to  have  been  ignited  and  a 
sustained  reaction  is  taking  place;  comparisons  between  numerical  and  experimental  data 
are  only  performed  for  cases  where  the  projectile  is  sufficiently  upstream  of  the  ignition 
point.  Simulation  of  mixture  ignition  will  be  included  in  future  reports.  The  global  reaction 
model  should  not  be  expected  to  reproduce  exact  profiles  of  chemical  species  nor  should  it 
necessarily  model  the  ignition  phase.  Rather,  the  results  obtained  with  the  global  model 
should  be  representative  of  the  overall  flowfield  with  regard  to  the  mixing  of  the  fuel  and 
the  overall  heat  release.  Therefore,  the  parameters  chosen  for  comparison  should  be  static 
temperature,  pressure,  distribution  of  mixing  and  heat  release. 

The  computational  domain  for  reacting  flow  simulations  extends  from  the  subsonic 
flow  upstream  boundary  (formed  by  the  normal  shock)  that  is  usually  positioned  on  the 
projectile  body  (Figure  3)  to  a  far-field  downstream  boundary  which  is  also  assumed  to  be 
subsonic.  For  the  thermally  choked  mode  the  downstream  boundary  eventually  becomes 
sonic;  however,  the  computational  domain  was  not  extended  to  this  point.  The  upstream 
flow  profiles  are  supplied  from  supersonic  non-reacting  flow  computations  (USA-PG  code) 
while  downstream  conditions  are  extrapolated.  Simulations  are  performed  in  the  projectile- 
fixed  reference  frame  (i.e.  moving  tube  wall)  and  subsequently  transformed  to  the  inertial 
frame.  Reacting  flow  simulations  are  performed  for  laminar  flow  modeling  (i.e.  fit  =  01 
while  results  for  reacting/turbulent  flow  will  be  included  in  future  reports. 

The  two  codes  -  USA-PG  and  BRL-LFD  RAMCOMB  -  can  be  used  in  either  an 
iterative  fashion  to  obtain  the  supersonic/subsonic,  inert/reacting  flowfield  or  in  a  single 
sweep  computation  where  the  flowfield  is  partioned  by  a  normal  shock  of  known  location. 
The  former  method  has  been  used  with  success  for  a  similar  mixed  flow  problem.11  In  the 
latter  method  the  supersonic/inert  flow  is  computed  using  the  USA-PG  code  which  supplies 
flow  profiles  for  subsonic/reacting  computations  using  RAMCOMB  after  the  normal  shock 
is  imposed.  As  an  initial  test  of  the  codes  the  single  sweep  method  was  used  for  results 
in  the  present  report  using  normal  shock  locations  derived  from  data  reported  by  the 
University  of  Washington.1 

The  computational  grid  for  the  38mm  configuration  2  starts  .1074m  (64%  body  length) 
from  the  projectile  nose  (i.e.  the  measured  normal  shock  position)  and  ends  .7m  down- 

33  Carnahan ,  B.,  Lather,  H.A.,  and  Wilkes,  J.O.,  Applied  Numerical  Methods,  John  Wiley  and  Sons,  Mew  )'ork,  1969. 
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stream.  In  the  streamwise  direction  201  points  are  distributed  and  51  points  are  used 
between  the  body  and  the  launch  tube  wall.  For  the  120mm  configuration  2  the  same  grid 
dimensions  are  distributed  from  .3441m  from  the  projectile  nose  (assumed  normal  shock 
location,  64%  body  length)  to  1.7m  downstream.  Calculations  were  performed  on  the  US 
Army  TACOM  and  US  Army  BRL  CRAY-2  supercomputers  in  about  .25  CPU  hours  per 
case. 


Figure  19  shows  computed  and  measured  launch  tube  wall  pressures  for  the  38mm 
configuration  2.  Measured  data  show  an  apparent  normal  shock  at  .1074m.  Flow  upstream 
of  the  normal  shock  was  computed  using  USA-PG  and  agreement  with  data  is  very  good. 
This  indicates  that  the  inert  and  reacting  flow  may  be  partitioned.  Downstream  flow  was 
computed  using  RAMCOMB  and  agreement  with  data  is  fair.  A  second  pressure  rise  after 
the  normal  shock  (x  ~  .2m)  is  not  reproduced  by  the  simulation  indicating  that  the  reaction 
zone  may  be  spread  over  a  larger  distance  than  represented  by  the  global  one-step  reaction 
with  Arrhenius  rate  terms  (see  comments  in  section  IV. 3).  In  addition,  the  measurements 
were  made  for  a  projectile  with  bore-riding  fins.  The  location  of  these  fins  with  respect  to 
the  launch  tube  wall  pressure  transducer  is  not  precisely  known  (i.e.  the  projectile  rotates 
in  the  tube).  Thus,  measurements  shown  in  Figure  19  may  represent  a  tube  wall/fin 
interaction.  Further  investigation  is  warranted.  The  computed  and  measured  pressures 
show  a  decreasing  trend  downstream  of  the  body  (x  >  .17m)  toward  a  sonic  condition, 
characteristic  of  subsonic  combustion. 

Figure  20  shows  a  comparison  between  computed  38mm  and  120mm  geometry  tube 
wall  pressures  for  reacting  flow.  Streamwise  distance  has  been  nondimensiomdized  by  the 
respective  body  lengths.  The  location  of  the  normal  shock  was  assumed  to  be  at  the 
same  location  (relative  to  total  body  length)  for  both  geometries.  Note  that  the  inert  flow 
pressures  (x/L  <  .64)  are  very  similar;  this  result  consistent  with  that  of  Figure  18.  The 
level  of  pressure  in  the  reaction  zone  is  smaller  for  the  120mm  geometry. 

Figure  21  shows  a  comparison  between  computed  38mm  and  120mm  geometry  body 
wall  pressures  for  reacting  flow.  No  measurements  have  been  reported  on  the  projectile 
body  to  date.  Pressure  oscillations  on  the  body  surface  (0  <  x  <  .15m)  are  caused  by 
pressure  signals  that  travel  upstream  in  the  near- wall  flowfield  and  may  be  related  to  the 
use  of  a  body-fixed  coordinate  system.  Pressures  are  similar  for  both  geometries  except 
for  a  pressure  rise  at  x/L  ~  .6  and  smaller  pressures  in  the  reaction  zone  for  the  120mm 
geometry.  Integration  of  these  pressures  yields  a  thrust  coefficient  (T / Poo-4)  of  1.4S  for  the 
38mm  geometry  and  1.26  for  the  120mm  geometry.  The  reference  quantities  are  tube  fill 
pressure  (Poo)  and  tube  cross-sectional  area  (A).  For  a  similar  fuel  mixture  in  the  3Smm 
system,  Reference  5  reports  a  thrust  coefficient  of  2.3.  In  dimensional  form  the  computed 
thrust  (T)  is  3400N  and  28700N  for  the  38mm  and  120mm  geometries,  respectively. 

Figure  22  shows  reaction  rate  contours  for  the  38mm  configuration  2.  The  dashed 
line  at  y  =  .019m  represents  the  launch  tube  wall  while  other  dashed  lines  represent  the 
boattail-base  body  geometry;  recall  that  the  calculation  is  started  at  x  =  .1074m  from  the 
projectile  nose.  Note  that  the  figure  is  plotted  with  a  skewed  scale  in  order  to  include  a 
larger  flowfield  area  around  the  body  thus,  body  angles  are  exaggerated.  The  clustering 
of  contour  lines  near  the  imposed  normal  shock  (x  =  .1074m)  as  well  as  tube  and  body 
walls  (i.e.  boundary  layers)  indicates  large  reaction  rate  gradients  at  these  locations.  The 
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magnitude  of  the  reaction  rate  is  also  largest  in  these  regions.  Flow  heating  in  these  regions 
is  the  primary  reason  for  increased  reaction.  Rate  gradients  are  reduced  in  the  projectile 
wake  where  reaction  products  accumulate  reducing  the  magnitude  of  the  reaction  rate  (see 
Equation  37)  as  shown  in  Figure  23.  These  simulations  indicate  that  a  larger  percentage  of 
the  reaction  (and  thus  high  pressure)  occurs  on  the  afterbody  as  opposed  to  the  projectile 
base.  Thus,  the  widely  used  schematic  shown  in  Figure  3  should  be  revised. 


V.  CONCLUSIONS  AND  FUTURE  WORK 


Computational  fluid  dynamics  solutions  of  the  Navier-Stokes  equations  have  been 
applied  to  both  non-reacting  and  reacting  in-bore  flowfields  for  a  thermally  choked  ram 
accelerator  projectile  launch  system.  Separate  computational  fluid  dynamics  codes  for  inert 
flow  (perfect  gas)  and  reacting  flow  (finite-rate  chemistry)  were  employed.  The  codes  were 
used  in  an  uncoupled  fashion  (i.e.  specified  normal  shock  location)  to  simulate  the  entire 
flowfield.  Good  comparison  between  computed  and  measured  pressures  for  non-reacting 
flow  in  a  38mm  system  were  achieved  including  the  investigation  of  viscous  modeling 
and  geometric  scaling  to  a  120mm  system.  Good  comparison  was  also  achieved  between 
computed  and  measured  pressures  in  the  subsonic  reaction  zone. 

Solutions  that  employ  a  coupling  technique  with  USA-PG  and  BRL-LFD  RAMCOMB 
are  in  progress  and  will  reported  at  a  future  date.  These  solutions  will  be  used  to  predict 
(and  compare  with  measured  data)  projectile  thrust  derived  from  a  variety  of  ram  accel¬ 
erator  fuel/oxidizer  mixtures  and  projectile  velocities.  Tests  conducted  at  the  University 
of  Washington5  have  shown  that  thrust  decreases  with  velocity,  reaching  a  minimum  near 
the  Chapman-Jouguet  detonation  speed,  and  then  showing  an  increasing  trend.  Numerical 
simulations  will  attempt  to  duplicate  this  trend.  Investigation  into  the  effects  of  turbu¬ 
lence  modeling  for  reacting  flows,  alternative  reaction  rate  expressions,  and  global  versus 
multi-step  reaction  mechanisms  is  warranted.  Parametric  studies  of  initial  gas  mixture  and 
pressure,  projectile  velocity  and  geometry  (38mm  and  120mm  scales)  and  other  effects  are 
also  planned. 

The  use  of  a  fully  coupled  single-code  CFD  capability  is  also  underway  at  the  BRL 
using  the  latest  versions  of  the  USA-series  codes.  These  codes  include  both  equilibrium 
chemistry  and  finite-rate  chemistry  fully  coupled  to  the  gasdynamics  and  arc  capable 
of  time-accurate  simulations.  The  USA-series  code  will  also  be  used  to  investigate  the 
oblique  detonation  and  superdetonative  stages  of  ram  accelerator  operation.  In  these 
stages  a  mixed-mode  flow  is  not  involved  since  combustion  occurs  primarily  in  supersonic 
flow.  Certain  aspects  of  numerical  simulations  will  therefore  be  greatly  simplified.  Three- 
dimensional  simulations  are  planned  and  will  include  the  bore-riding  fins  and  effects  of 
projectile  balloting  in  the  launch  tube. 
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Figure  1.  University  of  Washington  Ram  Accelerator  facility  and  typical  3$mm  projectile. 
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Figure  2.  Ram  Accelerator  in-bore  projectile  geometry  (excluding  bore-rider  fins);  lengths 
in  meters,  angles  in  degrees. 
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Figure  3. 


Schematic  of  Ram  Accelerator  in-bore  flowfield;  thermally  choked  mode. 
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Figure  9.  Computed  and  measured  tube  wall  pressures  for  non-reacting  flow.  3Smm 
Config.  2,  Moo  =  3.022,  =  1102m/s,  =  25atm. 
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Figure  10.  Computed  body  wall  pressures  for  non-reacting  flow.  3Smm  Confl 
A/oo  =  3.022,  Uoo  =  1102m/s,  =  25atm. 
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Figure  12.  Computed  tube  wall  pressures  for  non-reacting  turbulent  flow,  3Smm  Config. 
1  and  Config.  2,  M0 0  =  3.022,  U<x>  =  1102m/s,  =  25atm. 
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Figure  14.  Computed  tube  wall  pressures  for  non-reacting  flow,  120mm  Config. 
Moo  =  3.022,  Coo  =  1102m/s,  Poo  =  25atm. 
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Figure  15.  Computed  tube  wall  pressures  for  non-reacting  turbulent  flow,  120mm  Config 
1  and  Config.  2,  =  3.022,  =  1102m/s,  Px  =  25atm. 
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Figure  16.  Computed  tube  wall  pressures  for  non-reacting  turbulent  flow,  120mm  Config. 
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Figure  17.  Computed  tube  wall  pressures  for  non-reacting  turbulent  flow,  38mm  Config. 
2,  120mm  Config.  1  and  Config.  2,  M <*,  =  3.022,  Uoo  —  1102m/s,.Pco  =  25atm. 
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Figure  22.  Computed  reaction  rate  contours  for  reacting  flow,  38mm  Config.  2 
Moo  =  3.95,  Uoo  =  1430m/s,  =  20atm. 
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Figure  23.  Computed  product  mass  fraction 
Moo  =  3.95,  Uoo  =  1430m/s,  =  20atm. 
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List  of  Symbols 


A  cross  sectional  area  of  launch  tube 

cp  specific  heat  capacity,  constant  p 

cv  specific  heat  capacity,  constant  volume 

C  specific  reaction  rate  constant 

e  specific  total  internal  energy 

Ea  activation  energy 

F,G,H  flux  vectors  (Eq.  1) 
h  molar  specific  enthalpy 

h  total  enthalpy 

J  flow  rate  or  flux 

k  turbulence  kinetic  energy 

L  total  body  length 

m  species  mass  fraction 

M  Mach  number 

M.  molecular  weight  . 

n  stoichiometric  coefficient 

N  number  of  species 

p  static  pressure 

Pr  Prandtl  Number 

q  heat  transfer  rate 

r  radial  direction 

R  reaction  rate  per  unit  volume 

3?  specific  gas  constant,  (7  —  l)cp/7 

universal  gas  constant,  ?  M j 

Re  Reynolds  Number 

Sc  Schmidt  Number 

t  time 

T  static  temperature 

T  thrust 

u  axial  velocity  (USA-PG  code);  radial  velocity  (RAMCOMB  code) 

U  mean  streamwise  velocity 

v  radial  velocity  (USA-PG  code);  azimuthal  velocity  (RAMCOMB  code) 

V  magnitude  of  the  local  velocity  vector 

V  uf  +  v6  +  wz 

w  velocity  component  in  the  axial  direction 

W  dependent  variable  vector  (Eq.  1) 

x,  y  cartesian  coordinates 

z  axial  direction 

A Hj  heat  of  formation 

A  Hr  heat  of  reaction 

Greek  Symbols 

7  ratio  of  specific  heats,  cp/ c„ 
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r  diffusion  coefficient 

8  stoichiometric  air/fuel  mass  ratio 

t  turbulence  dissipation  rate 

t /  transformed  coordinate 

9  azimuthal  direction 

H  molecular  viscosity 

£  transformed  coordinate 

p  density 

0j 

0  k  /Xeff/rfc 

a xx,  <yTr  normal  stress  tensors 
r  transformed  time 

txt  shear  stress  tensor 

f  shear  stress  vector 

(f>  general  flow  variable 

rp  streajn  function 

u  vorticity 

Superscripts 


unit  vector 

total  or  stagnation 

rate 


Subscripts 


eff  effective 

/  adiabatic  flame  quantity 

h  enthalpy 

1  initial  or  base-state  quantity 

j  j- th  mixture  component  or  species 

k  turbulence  kinetic  energy 

p  constant  pressure 

r  radial  component  or  radial  direction 

t  turbulence  quantity 

v  constant  volume 

x  x-direction 

2  axial  component 

e  turbulence  dissipation  rate 

r]  ^-direction  transform  coefficient 

9  azimuthal  component 

£  ^-direction  transform  coefficient 

oo  freestream  quantity 
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